#rm(list=ls()) 

library(SDMTools) ; library(maptools) #load the necessary libraries

#to make target group background, have to read in all of the species files and get all the unique occurrences for all the species in that class
	spp.dir = "/home/ctbccr/Ramona/Species_data/" ; setwd(spp.dir)
	climate.dir = '/home/jc165798/Climate/CIAS/Australia/1km/baseline.76to05/bioclim/' #define the climate directory
	#climate.dir = "/home/ctbccr/Ramona/Climate_Data/bioclim_mxe/"
	taxa = c("Amphibians", "Birds","Mammals","Reptiles","Plants") ; tax= taxa[1]
	
	occur = read.csv("myamphiok.csv") ; occur = occur[,c(1:3)] #read in the observations 
	#creating background base file for every species
	base.asc = read.asc.gz("/home/jc165798/Climate/CIAS/Australia/1km/bioclim_asc/RCP3PD_cccma-cgcm31_2015/bioclim_01.asc.gz")
	pos = as.data.frame(which(is.finite(base.asc),arr.ind=TRUE))
	pos$lat = getXYcoords(base.asc)$y[pos$col] #extract the longitudes
	pos$lon = getXYcoords(base.asc)$x[pos$row] #extract the latitudes
	
	#for future reference: if it says output/maxentResults.csv doesn't exist - it might be because the long/lat files are the wrong way around!!
	clim.vars = paste('bioclim_',sprintf('%02i',c(1,4,12,16,17)),sep='') #define the climate variables of interest
		for (clim.var in clim.vars) { cat(clim.var,'\n') #cycle through each of hte climate variables
		occur[,clim.var] = extract.data(cbind(occur$longitude,occur$latitude), read.asc(paste(climate.dir,clim.var,'.asc',sep=''))) #append the climate data
		#repeat for substrate variables
		}
	occur = na.omit(occur[,c('species','longitude','latitude',colnames(occur)[grep('bioclim',colnames(occur))])]) #reorder the occur
	
	env.vars = c("clay.new", "hpedality.new","ksat0.new")
	env.dir = "/home/ctbccr/Ramona/Substropo_data/new/"
	#extracting substrate variable for occurrence file
		for (env.var in env.vars) { cat(env.var,'\n')
			occur[,env.var] = extract.data(cbind(occur$longitude, occur$latitude), read.asc(paste(env.dir,env.var,'.asc',sep='')))
			#repeat for substrate variables
			}
	#define the background
	#write.csv(bkgd,'all_Aust_bkgd.csv',row.names=FALSE); rm(bkgd) #write out the file and remove it from memory
	#now we repeat the extraction for the background file
	colnames(pos) = c("row","col","latitude","longitude") ; bkgd.all = pos ; bkgd.all = bkgd.all[c(1,2,4,3)]
	clim.vars = paste('bioclim_',sprintf('%02i',c(1,4,12,16,17)),sep='') #define the climate variables of interest
		for (clim.var in clim.vars) { cat(clim.var,'\n') #cycle through each of hte climate variables
		bkgd.all[,clim.var] = extract.data(cbind(bkgd.all$longitude,bkgd.all$latitude), read.asc(paste(climate.dir,clim.var,'.asc',sep=''))) #append the climate data
		#repeat for substrate variables
		}
	#bkgd.all = na.omit(bkgd.all[,c('species','longitude','latitude',colnames(bkgd.all)[grep('bioclim',colnames(bkgd.all))])]) #reorder the occur
	#extracting substrate variable for occurrence file
		for (env.var in env.vars) { cat(env.var,'\n')
			bkgd.all[,env.var] = extract.data(cbind(bkgd.all$longitude, bkgd.all$latitude), read.asc(paste(env.dir,env.var,'.asc',sep='')))
			#repeat for substrate variables
			}
	bkgd.all$species = 'bkgd'	#; bkgd.all = unique(bkgd.all)
	bkgd.all = bkgd.all[c(ncol(bkgd.all),4:(ncol(bkgd.all)-1))] ; bkgd.all = na.omit(bkgd.all)
	
	write.csv(bkgd.all,"bkgd.all.csv",row.names=F)
	
	species = unique(occur$species)
	taxa.folder = paste('/home/ctbccr/Ramona/SDM/',tax,"/",sep='') ; wd = taxa.folder  ; setwd(wd)
	
	for (spp in species[1:15]) {	
	#spp=species[2]	
		spp.dir = paste('/home/ctbccr/Ramona/SDM/',tax,"/models/",spp,'/',sep='') ; dir.create(spp.dir)#define the species directory
		#if (!file.exists(paste(spp.dir,'output',sep=''))){
		dir.create(paste(spp.dir,'output',sep=''),recursive=TRUE) #create the species directory 
		toccur = occur[occur$species==spp,] #; colnames(occur) = c("SPPCODE","Lon","Lat","bioclim_01", "bioclim_04", "bioclim_05","bioclim_06", "bioclim_12","bioclim_15","bioclim_16","bioclim_17" )
		#coordinates(toccur) = ~Lon+Lat #get a subset of the data
		write.csv(toccur,paste(spp.dir,'occur.csv',sep=''),row.names=FALSE) #write out the file and remove it from memory

		#to make a background limited to subbioregion	
	# = read.csv('bkgd.au.csv') # define directory
	clipasc = read.asc('/home/ctbccr/Ramona/IBRAs/ibra7.asc')
	regions = extract.data(toccur[,2:3],clipasc) #extract the regions occurrences fall within
	regions = unique(regions) #identify the unique regions
	posoi = which(clipasc %in% regions) #identify locations within these regions
	clipasc[which(is.finite(clipasc))] = 0; clipasc[posoi] = 1 #change everything to 0 and only locations in bioregions of interest to 1
	write.asc.gz(clipasc,paste(spp.dir,spp,'.bioregion.clip.asc',sep='')) #write out the clipping ascii
	#pos = pos[c(1:4)]
	bkgd.all$clip = extract.data(cbind(bkgd.all$longitude,bkgd.all$latitude),clipasc) #extract the clip information
	pos.of.interest = bkgd.all[which(bkgd.all$clip == 1),] ; 
	trows=sample(nrow(pos.of.interest),size=10000,replace=T)
	bkgd = pos.of.interest[trows,]
	bkgd.spp = bkgd[,-ncol(bkgd)] ; bkgd.spp = na.omit(bkgd.spp) #; bkgd.spp = bkgd.spp[,c(ncol,3:(ncol
	write.csv(bkgd.spp,paste(spp.dir,spp,"_bkgd.csv",sep=''),row.names=FALSE)
	write.csv(bkgd.spp,paste(spp.dir,"bkgd.csv",sep=''),row.names=FALSE)
	
		#maxentResults.csv = paste(spp.dir,"output/maxentResults.csv",sep = '') ; cv = paste(spp.dir,"output/maxentResults.crossvalide.csv",sep = '') ; maxentResults.csv = cv = NULL
		zz = file(paste(spp.dir,'01.create.model.sh',sep=''),'w') #create the shell script to run the maxent model
		cat('#!/bin/bash\n',file=zz)
		cat('cd $PBS_O_WORKDIR\n',file=zz)
		cat('source /etc/profile.d/modules.sh\n',file=zz)
		cat('module load java\n',file=zz)					
		cat('java -mx1024m -jar ',wd,'maxent.jar -e ',spp.dir,'bkgd.csv -s ', spp.dir, 'occur.csv -o output nothreshold nolinear noquadratic noproduct noautofeature nowarnings novisible replicates=10 nooutputgrids -r -a \n',sep="",file=zz)
		cat('cp -af ',spp.dir,',output/maxentResults.csv ',spp.dir,',output/maxentResults.crossvalide.csv\n',file=zz)
		cat('java -mx1024m -jar ',wd,'maxent.jar -e ',spp.dir,'bkgd.csv -s ', spp.dir, 'occur.csv -o output nothreshold nolinear noquadratic noproduct noautofeature nowarnings novisible nowriteclampgrid nowritemess writeplotdata -P -J -r -a \n',sep="",file=zz)
		close(zz)
		setwd(spp.dir); system(paste('qsub -m n -N ',spp,' 01.create.model.sh',sep='')); setwd(wd) #submit the script
		}
	
#nolinear noquadratic noproduct noautofeature
